Self-Emerging and Turbulent Chimeras in Oscillator Chains 
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We report on a self-emerging chimera state in a homogeneous chain of nonlocally and nonlinearly 
coupled oscillators. This chimera, i.e., a state with coexisting regions of complete and partial 
synchrony, emerges via a supercritical bifurcation from a homogeneous state and thus does not 
require preparation of special initial conditions. We develop a theory of chimera based on the 
equations for the local complex order parameter in the Ott-Antonsen approximation. Applying a 
numerical linear stability analysis, we also describe the instability of the chimera and transition to 
phase turbulence with persistent patches of synchrony. 
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Populations of coupled oscillators is a paradigmatic 
model of nonlinear science, with numerous applications 
from purely physical ones like Josephson junction arrays 
and coupled lasers to biologically and even socially im- 
portant pp. While in globally coupled ensembles and in 
networks one is mostly interested in the features of syn- 
chronization and desynchronization, spatially extended 
oscillating systems demonstrate a variety of pattern- 
forming phenomena. One of the most spectacular recent 
findings are the so-called chimera states (CSs) which are 
observed in otherwise completely synchronizable oscilla- 
tory media if the system starts from specially prepared 
initial state. CSs are characterized by the coexistence 
of regions with locally synchronized oscillators and re- 
gions where the oscillators phases are not locked but yet 
not completely incoherent. CSs were initially discovered 
and explained theoretically in [2] , and then received more 
analytical treatment in [3]. Following those pioneering 
works on CSs, a large body of observations and analy- 
sis of similar regimes has been recently published, see [4] 
and references therein. 

In this Letter, we add another species to the zoo of 
chimeras. The crucial difference is that our CS does not 
require special initial conditions: It emerges from a gen- 
eral initial state and is thus denoted as self-emerging. 
CS is stable close to the bifurcation, but with a fur- 
ther variation of the parameter it becomes turbulent, so 
that synchronous and partially synchronous patches in- 
termingle irregularly. The key elements of our model are 
Stuart-Landau oscillators, coupled through an exponen- 
tially decaying kernel as in the original chimera setup, 
but with a difference that the coupling is nonlinear in 
the sense of [5]. First, we numerically demonstrate the 
existence of CS and then explain it in the phase dynam- 
ics framework with the help of reduced equations for the 
local order parameter. Our main theoretical tools are 
the equations for the complex order parameter in the so- 
called Ott-Antonsen (OA) approximation [6|. It exploits 
a parametrization of the probability density for ensem- 



bles of sinusoidally coupled phase oscillators and results 
in a closed equation for the order parameter. The OA 
ansatz is closely related to the Watanabe-Strogatz the- 
ory [7] which is exact but does not yield closed equations 
in terms of the order parameter. A connection between 
these two theories has been established in [8]. 

As a basic model we consider a one-dimensional, peri- 
odic in space chain of the length L = 2£ of nonlocally 
coupled Stuart-Landau oscillators 



d t A= (l + iu)A-\A\ 2 A + eZ , 



(1) 



where A = A(x, t) is the complex amplitude, oj is the 
natural frequency of the oscillators, Z = Z (x,t) is the 
coupling force acting on the oscillator at and e is a 
small coupling constant. The coupling is organized via a 
convolution of A (x, t) with the weight function G (x) = 
ce-N : 



B(x,t) 



dx'G(x-x')A(x',t) , 



(2) 



where the constant c ensures that J^ £ G (x) dx = 1. The 
forcing is then defined as 



Z = e [ ^ B »B. 



(3) 



The phase shift (3 (\B\) = f3 + fii\B\ 2 , chosen in spirit 
of p)], accounts for possible nonlinearity effects in the 
coupling, i.e., the forcing Z in Eq. depends on the 
higher-order powers of B. The weight function G (x) is 
kept fixed, so that by variation of L we change the ratio 
between the coupling wi dth and the system size. 

We integrated Eqs. (IJ2) using 2 17 sites in x with 
e = 0.01, cD = 0, A) = 0.4tt and = (tt/2 - /3 ) /0.36 
(this choice of j3\ will be apparent later on) using the 
Runge-Kutta 4th order scheme. Initial conditions were 
chosen close to the completely desynchronized state. To 
measure the synchronization between neighboring oscil- 
lators on a mesoscopic scale, we calculated the coarse- 
grained A(x,t) by averaging A(x,t) over 2 10 closest 
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FIG. 1. (color online) Results of simulation of Eq. |T]). (a) and 
(b) show stationary patterns, (c-e) show a turbulent one. In 
(a-c) we show profiles of \A\ (dotted red line), \A\ (solid green 
line), and \B\ (long-dashed blue line). The system size L is 
specified in each frame. In (b), thin pink solid line (largely 
occulted by the green one) shows a snapshot from simula- 
tion with specially prepared initial state, see text. In panels 
(d) and (e), showing space-time plots of \A\ and \B\, respec- 
tively, red color (darker in b/w print) encodes values close to 
1, blueish color (brighter in b/w print) encodes values close 
to 0. In (d) and (e), the time span of simulations is 5000 
dimensionless time units and L = 7.5. 



neighbors. We varied the system size L and looked at 
the dynamics of \B\ and \A\. For small L, we observed a 
spatially homogeneous, uniformly rotating, and partially 
synchronous state with \A\ = \B\ ~ 0.6, like in globally 
coupled ensembles [5 j. This state becomes unstable if the 
system size L exceedes the critical value L c w 5.1. For 
L > L c , a spatially modulated profile of \B\ emerged, see 
Fig. [I] (a). The profile of \A\ was also stationary up to 
finite-size fluctuations. Close to the transition, synchro- 
nization was only partial, with \A\ < 1 for all x. Howevei\ 
for L > 5.35, the profile of \A\ reached unity (see Fig. Ill 
(b)), which means that locally all oscillators in that re- 
gion are completely synchronized, whereas in the regions 
with \A\ < 1 the local synchronization was only partial. 
Such state is a stationary chimera [2j [3]. With further 
increase of L, this regime becomes unstable and evolves 
into a turbulent state where synchronized patches with 
\A\ « 1 appeared at random places and disappeared af- 
ter some lifetime, see Fig. [I] (c)-(e). We call this state 
a turbulent chimera. Below, we present a theoretical de- 
scription of stationary CSs in terms of a reduced phase 
model. 

For a small coupling constant e we can assume that 
the amplitude of each Stuart-Landau oscillator remains 
unperturbed and the system can be described solely by 
its phase 6(x,t). In the thermodynamic limit, we de- 
scribe the population with a probability density distri- 
bution p = p(x,0,t). This density obeys the continu- 



ity equation d t p + &q (pv) = with the velocity field 
v (x, 0, t) = = uo + Im (Z (x, t) e-' l °) , where Z (x, t) is 
defined in Eq. |3]). (We rescaled the coupling parameter 
e — >• 1 by rescaling time and frequency t — »• et, uo — uo/s). 
Next, we define the complex order parameter by 

z(x,i) = re i(p = / d0e w p (x, 6, t) , 
Jo 

so that B (x,t) is expressed as 



B(x,i) = b(x,t)e 



f dx'G(\ 
J-t 



x — x'\) z (x f , t) . 

(4) 

Now we apply the theories developed in [6]-[8] and use the 
OA equation for the complex order parameter z(x,i) 



d t z 



(Z -z 2 Z*) 



(5) 



Together with Eqs. (3j4) it constitutes a closed system 
which we analyze to study pattern formation in the do- 
mains of various length L. This system is similar to 
Eqs. (13,14) in [4] with the exception of the nonlinear 
phase shift /?. 

The applicability of the OA theory requires justifica- 
tion, since for identical oscillators under a common forc- 
ing the ensemble dynamics can go beyond the OA the- 
ory [8]. For non-identical oscillators, the phase distri- 
bution function asymptotically fulfills the conditions for 
the OA theory (one says that a solution is attracted to 
the invariant OA manifold) [6]. Here, even though the 
oscillators are identical, the force that acts on them is in- 
homogeneous in space (see profiles of \B\ in Fig.[T]), which 
plays the same mixing role as the non-identity of oscilla- 
tors in [6 and the phase distribution will tend to the OA 
manifold as well. To support this claim numerically, we 
ran an additional simulation with initial condition that 
strongly violated the OA ansatz. This simulation also 
resulted in a stationary CS. The profile of \B\ was in- 
distinguishable from the previous simulation with nearly 
uniform initial distribution of phases. The profile of \A\ 
(thin solid pink line in Fig. [l] (b)) is close to the previ- 
ous one up to the finite-size fluctuations in the region of 
small | B | values. The differences in \A\ are mostly pro- 
nounced in the domain with nearly uniform forcing \B\, 
what agrees with the aforementioned reasoning. We also 
followed the deviation from the OA manifold looking at 
the quantity |(e l6> ) 2 — (e l26> )| ((•) denote averaging over 
small x- neighbor hood as described above), which van- 
ishes on the OA manifold. We found that even for initial 
conditions off the OA manifold, asymptotically it always 
descreased in time. 

With a proper choice of /3q and Eq. (|5| supports 
three types of homogeneous steady states: (i) The com- 
pletely desynchronized one zo — 0, (ii) the fully synchro- 
nized one z\ = e lQlt with Qi = uj + sin (/3q + and 
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FIG. 2. The real part of the stability eigenvalues A of z soq _ 
in dependence on k (the wavenumber of perturbation). k c ~ 
27r/5.10 denotes the critical wavenumber. 



(iii) the intermediate regime of self-organized quasiperi- 



odicity (SOQ) [5J 



with r s , 



and fL 



■Hi- 



1 ' soq/ 

and z S( 



' soq c vviuii / SO q — A/ (3 1 

Due to the normalization of 
the kernel G, zo, z \-> an d ^soq persist for all £. Knowing 
from [5] that both zo and zi are unstable, we focus on 
z soq . A linear stability analysis of z soq with the pertur- 
bation oc y e xt + lkx results in the eigenvalue problem 



ye 
y*e 



ikx 
ikx 



^12 



Cl2 
^11 



ye 
y*e 



ikx 



where Cn = (i (r 2 oq - l) + a (fc, £) (i 
and C12 = a (k,£) (irf oq - fa (l - 
dependent wavenumber factor 



• 2 )) 

soq// 



(i 



,2 

soq V" 1 " ' soq 

with the 



))) 



a(k,£) 



1 



(7c sin (H) -cos(H)) 



(1 



0(l + fc 2 ) 



(6) 



which converges to (l + k 2 ) as £ — » oo. The eigenval- 
ues Ai 5 2 of matrix £ depend on both fc and On the 
infinite domain ^ = oo, we find that ReAi (fc, oo) > 
and Im Ai 5 2 oo) = for < k < k c with the crit- 
ical wavenumber k c = \ R I^ 2/3 ° /0 • With a finite £. 
the critical wavenumber can differ from k c since £ en- 
ters a(k,£). The critical system length L c = 2£ c at 
which the instability occurs, is determined by the condi- 
tion Re A (tt/£ c , £ c ) = 0. To exemplify the instability, we 
chose f3 = 0.47T and f3\ such that /3 + Pi x 0.6 2 = 7r/2, 
so that r SO q — 0.6 (cf. numerics of Fig. [I}. Solving 
the eigenvalue problem, we obtained the critical length 
L c « 5.09 in a nice correspondence with the results of 
the direct numerical simulations. 

We expect that the instability described above results 
in uniformly rotating spatially inhomogeneous regimes, 
thus we look for a solution z c h (x, t) = r (x) e^+^fa)] 
with some unknown frequency Q and spatial functions 
r (x) and ip(x). For the forcing we assume the ansatz 
B(x,t) = b (x) e i[nt+*P(x)) _ when substituting this in 
Eq. (§, we have to differentiate two cases: (i) for re- 
gions where r < 1 the stationarity/stability conditions 
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FIG. 3. (color online) Stationary profiles of r for different sys- 
tem sizes L. On the left, r = \z\ (solid red line) and b (dashed 
blue line) in comparison to A (horizontal short-dashed line) 
are depicted, the corresponding stability eigenvalues of the 
Jacobian matrix are on the right. 



yield 

r (x) 
ip(x) 

with the detuning A 
we have 



x/>(x)+ P(b(x))- ir/2, 



(7) 



£1 — cj; (ii) for regions with r — 1 



(x) = ijj (x) + /3 (b (x)) + arcsin T (x) . 



(8) 



where in both Eqs. ( pf|8] ) T (x) = A/b(x). Equation 
or ([§ plus Eq. @ constitute a closed system, which we 
could not solve analytically. Instead, we found the sta- 
tionary CSs in the frame rotating with frequency em- 
ploying the Newton's method for Eq. ^ discretized in 
x. In the Newton's calculations, an additional phase pin- 
ning condition was imposed in order to pick up a unique 
solution from the family of phase rotations. The number 
of unknowns and equations was then balanced by taking 
Q as an additional unknown. After the solution z c h (x) 
had been found with a desired numerical accuracy, we 
looked at its linear stability by computing numerically 
the eigenvalues A of the Jacobian matrix evaluated at 
z c h(x). Eigenvalues with positive real parts would sig- 
nalize an instability of z c h (x). 

For smaller L close to L c , the profile of z c h (x) shows a 
moderate sine-like modulation, see Fig. [3] (a). For lar ger 
L, the profile of z c h (x) touches unity: A patch of com- 
plete synchronization emerges; this is a genuine CS. This 
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FIG. 4. (color online) Space-time plots of the coarse-grained 
order parameter \z\ (a) and of b (b) obtained by numerical 
integration of Eq. (|5| with L = 7.5. Red color (darker in b/w 
print) represents larger values close to one, bluish (brighter in 
b/w print) represents smaller values close to zero. We show 
simulation results from 1500 to 3000 dimensionless time units. 
The discretization in space x was done on 1024 sites, (c) Final 
snapshot of integration from (a) and (b): Solid red line shows 
non-coarse-grained r = \z\, green short-dashed line shows the 
magnitude of the coarse-grained z, and long-dashed blue line 
represents profile of b. 



occurs in the region with b > A, see Fig. [3] (c). For our 
choice of /?o and /3i, the regimes with r < 1 for all x are 
linearly stable (see Fig. [3] (a) and (b)), whereas the CS 
with a synchronized patch with r = 1 becomes unsta- 
ble as L increases, cf. Fig.[3](c-f). We believe that this 
instability is inherited from the instability of z±, due to 
the presence a plateau of complete synchronization with 

1*1 = 1. 

Our numerical simulations show that beyond the in- 
stability of the stationary CS, a turbulent state occurs, 
where a synchronous patch is persistent but appears at 
different places. We call this state a turbulent chimera. 
Figures [4] (a) and (b) show a result of numerical integra- 
tion of our system for L = 7.5. Since the order parameter 
z should be a coarse-grained quantity, in numerical simu- 
lations one has to average the numerically obtained field 
over a small spatial interval (in our simulation over 16 
neighboring sites) to get z. In Fig. [4] (a) the red patches 
with coarse-grained \z\ close to one show a larger de- 
gree of local coherence of oscillators in comparison to the 
blueish rest. These coherent regions persist, allowing us 
to characterize the irregular state as a chimera turbu- 
lence. The forcing; mass; nitude b (cf. Fig. [3] (b)) does 
not show any saturation, irregularly oscillating between 
smaller (blue color) and larger (red color) values. A typ- 
ical final snapshot of the system is shown in Fig. [3] (c). 
There is a fully synchronous patch near x/L « 0.5 where 

The presented results can be extended to non-identical 
oscillators with a Lorentzian distribution of frequencies. 



In this case the OA ansatz yields Eq. ([5| with an addi- 
tional damping term (so that the eigenvalue spectra of 
CS move from the imaginary axis into the left complex 
half-plane), which results in very similar states with the 
difference that synchronization is never complete. 

Summarizing, we have demonstrated the existence of 
chimera-like solutions in chains of nonlocally and nonlin- 
early coupled oscillators. In our system, chimera appears 
as a result of a long-wave instability via a supercritical 
bifurcation. Unlike the previously known CSs, we do not 
need to prepare initial conditions carefully to avoid com- 
plete synchrony, because it is unstable. We demonstrated 
that CSs are asymptotically well described within the 
OA theory even if the initial condition is not on the OA 
manifold. As the stationary chimera becomes unstable, 
it evolves into a turbulent one. In this novel regime the 
complex order parameter changes irregularly in space and 
time, but nevertheless the patches of synchronized oscil- 
lators appear persistently. This characterizes the chimera 
as an important pattern in the dynamics of nonlinearly 
coupled oscillators. 
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